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ABSTRACT 

The magnetorotational instability (MRI) is believed to be an efficient way to trans- 
port angular momentum in accretion discs. It has also been suggested as a way to 
amplify magnetic fields in discs, the instability acting as a nonlinear dynamo. Recent 
numerical work has shown that a large-scale magnetic field, which is predominantly 
azimuthal, axisymmetric and has zero net flux, can be sustained by motions driven by 
the MRI of this same field. Following this idea, we present an analytical calculation 
of the MRI in the presence of an azimuthal field with a non-trivial vertical structure. 
In the limit of small vertical wavelengths, we show that magnetorotational shearing 
waves have the form of vertically localized wavepackets that follow the classical MRI 
dispersion relation to a first approximation. We determine analytically the spatiotem- 
poral evolution of these wavepackets and calculate the associated mean electromotive 
force (EMF), which results from the correlation of the velocity and magnetic field per- 
turbations. The vertical structure of the azimuthal field results in a radial EMF that 
tends to reduce the magnetic energy, acting like a turbulent resistivity by mixing the 
non-uniform azimuthal field. Meanwhile, the azimuthal EMF generates a radial field 
that, in combination with the Keplerian shear, tends to amplify the azimuthal field 
and can therefore assist in the dynamo process. This effect, however, is reversed for 
sufficiently strong azimuthal fields, naturally leading to a saturation of the dynamo 
and possibly to a cyclic behaviour of the magnetic field. We compare these findings 
with numerical solutions of the linearized equations in various approximations, and 
show them to be compatible with recent nonlinear simulations of an MRI dynamo. 
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INTRODUCTION 



The magnetorotational instability (MRI) is believed to be r esponsible for turbulent motio n and angular momentum transport 
in accretion discs, at least those that are sufficiently ionized ( Balbus fc Hawlevlll99ll . 1998|). In its sim plest version it appears as 
a linear instability of a rotating shear flow in the presence of a uniform magnetic field |Velikhovlll959h . The nonlinear outcome, 



at su fficiently large Reynolds and magnetic Reynolds numbers, is magnetohydrodynamic (MHD) turbulence (jHawlev et al 



19951 ). More importantly, the MRI is also believed to act as a dynamo in accretion discs, meaning that the inductive effect 
of the motions driven by the instability is able to sustain the magnetic field aga inst resistive decay and so to perpetuate 
the turbulent motion even in the absence of an externally imposed magnetic field l Hawlev et al. 19961 ). The dynamo is fully 
nonlinear because the Lorentz force is essential to the operation of the MRI. 



The MRI dynamo can be studied in the well known shearing-box model (lHawlev et al.lll995l ). which is a local represen- 
tation of a differentially rotating disc, with boundary conditions that conserve the box-averaged magnetic field. When the 
initial conditions include a magnetic field with zero volume average, the MRI can be initiated but it has the opportunity to 
annihilate the magnetic field and thereby suppress the turbulence. In order for the magnetic field and the turbulence to be 
sustained indefinitely, dynamo action must occur. 

Recent numerical simulations of the MRI with zero box-averaged magnetic field have shown that the outcome depends 
in an important way on the numerical method and resolution unless explicit diffusion coefficients (viscosity and resistivity) 
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are included and the relevant dissipative scales are resolved in the calculation (jFromang fc Papaloizoul 120071 ; iFromang et al 



2007). These results cast doubt on the operation and efficiency of the nonlinear dynamo in accretion discs, especially in the 
regime of small magnetic Prandtl number (Pm). It is possible that the shearing box is too restrictive a model to describe 
the behaviour of accretion discs, as it conserves the box-averaged field. Different vertical bound ary conditions that allow 
horizontal flux to escape may assist in the operation of the MRI dynamo (jBrandenburg et al,lll995h . However, in the light of 
the results described below, we think this model is a good way to isolate the dynamo process from other effects, and can be 
seen as a minimal setup in which to study this process. 

Steady solutions of the M HD equations repres enting a self-sustaining MRI dynamo process have been obtained by nu- 
merical continuation methods IjRincon et al.l 120071) in a rotating shear flow between conducting walls. The magnetic field 
has zero net flux and is predominantly toroidal (azimuthal). In a condition of marginal stability, the MRI generates steady 
nonaxisymmetric motion. This induces a poloidal (radial and vertical) magnetic field which, through the action of Keplerian 
shear, sustains the toroidal field. These solutions reveal the workings of the dynamo but so far they have been found only for 
large Pm, perhaps because the method is restricted to steady states. 

We have recently identified a nd analysed a cyclical behaviour of the large-scale magnetic field in an MRI simulation in 
a shearing box with zero net flux ( Lesur fc O gilvie] 200S ). Th e processes involved in the cycle are related to those occurring 



in the wall-bounded steady solutions of 



Rincon et al 



(2007) and are probably relevant to the operation of the nonlinear 
dynamo in accretion discs. However, there is an importa nt role for shearing waves, whi ch are the typical nonaxisymmetric 
solutions obtained in the shearing box without any walls ( Goldreich fc Lvnden-Bell 19651 ). The nonaxisymmetric instabilities 
are transient in nature and difficult to analyse, yet they appear to be responsible for sustaining the magnetic field and limiting 
its growth. 

A common feature of these calculations is that the magnetic field is dominated by a large-scale, axisymmetric, azimuthal 
component. We can think of this as a 'mean' field if we consider a process of azimuthal averaging. At a simple level the dynamo 
operates because the MRI generates radial field from azimuthal, while shear creates azimuthal field from radial. However, 
the first process requires detailed understanding. In linear theory, the MRI of an axisymmetric azimuthal field gives rise to a 
nonaxisymmetric radial magnetic field perturbation that averages to zero because of its wavelike form. Since we are looking 
for an axisymmetric effect, we are interested in the average second-order effect resulting from the correlations of the velocity 
and magnetic field perturbations. This is described by the mean electromot ive force (EMF) of the MRI solutions. 

We interpreted the results of the simulations in Lesur fc Ogilvie (2008) by carrying out linear numerical calculations of 
the evolution of shearing waves in the presence of a large-scale azimuthal field with a sinusoidal dependence on the vertical 
coordinate, B x oc cosfcz. The waves undergo transient growth as a result of the MRI. We analysed the mean EMF associated 
with the waves and its feedback on the large-scale field. The radial EMF was found to have a resistive effect, reducing the 
energy of the azimuthal field. The azimuthal EMF £ x was found to generate a radial field which, in combination with the 
Keplerian shear, also affects B x . When B x is relatively weak, £ x results in an amplification of the field; when it is stronger, 
the sign of £ x is reversed and B x is reduced by the process. This leads to a cyclical behaviour in which the azimuthal field 
strength oscillates irregularly about a characteristic value. 

The purpose of this paper is to make an analytical and more general investigation into the behaviour of shearing waves 
in the presence of a vertically non-uniform azimuthal magnetic field. We aim to explain the systematic behaviour of the 
EMF, which is an essential part of the dynamo process; it is not our objective here to give a complete explanation of the 
dynamo. In order to make analytical progress we examine an asymptotic regime in which a separation occurs between the 
vertical length-scales of the large-scale field and of the shearing wave. We develop a general theory of vertically localized 
waves in Section 2 and derive expressions for the associated EMF. In Section 3 we compare these with the results of numerical 
calculations of the linearized equations. We summarize our findings and draw conclusions in Section 4. 



2 VERTICALLY LOCALIZED SHEARING WAVES 
2.1 Basic equations 

As a local description of an accretion disc, we adopt the model of the shearing sheet. We therefore consider a Cartesian 
coordinate system rotating with angular velocity Q e z , in which the basic state consists of a linear shear flow U — Sy e x . The 
radial, azimuthal and vertical directions therefore correspond to y, —x and z respectively, while the shear rate in a Keplerian 
disc would be S = (3/2)0. 

The effects we intend to describe do not depend on compressibility, stratification or diffusive processes. We therefore work 
with the ideal magnetohydrodynamic (MHD) equations for a homogeneous incompressible fluid. For simplicity, we also adopt 
units such that p = po = 1, which means that the magnetic field is expressed as the equivalent Alfven velocity. 

In addition to the linear shear flow, we include in the basic state a non- uniform azimuthal magnetic field B = B(z) e x of 
arbitrary form, which is responsible for the magnetorotational instability (MRI). The balance of forces in the basic state can 
be maintained by an appropriate pressure gradient. 
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We next consider small disturbances in the form of shearing waves in which the velocity and magnetic perturbations are 

Re [u(z, t) exp(ik x x + ik y y)] , Re [b(z, t) exp(ik x x + ik v y)] , (1) 
where k x = constant is the azimuthal wavenumber and 

ky — kx S ij 

is a time- dependent radial wavenumber that passes through zero at the chosen origin of time. The wave is leading for t < 
and trailing for t > 0. The linearized MHD equations are then 

ii x = (2fi — S)u v — ikxp + ik x Bb x + b z d z B, (3) 

ii y = —2Q,u x — ikyp + ik x Bb y , (4) 

iiz = — d z p + ik x Bb z , (5) 

b x = Sb y + \k x Bu x — u z d z B, (6) 

by — \k x Bu y , (7) 

b z = \k x Bu z , (8) 

ik x u x + ikyUy + d z u z = 0, (9) 

lk x b x + ikyb y + d z b z = 0, (10) 

where p represents the total pressure perturbation in the shearing wave. In fact, the last equation is required only as an initial 
condition, because its time-derivative is implied by the preceding equations. 

From the above we can obtain an equation for the azimuthally averaged energy density of the wave, 

i<9 t (i|u| 2 + i|6| 2 ) + \d z Ke{u z p) = iSRe(b* x b y - u x u y ) + \{d z B) Re{u x b z - u z b x ). (11) 

The factor of s here derives from the complex notation of equation ([1]), in which the average of the square of the velocity 
perturbation is ^|u| 2 , etc. When the energy equation is vertically integrated with appropriate boundary conditions, the term 
involving the pressure perturbation is eliminated. Potential sources of energy for the wave appear on the right-hand side of 
the energy equation. These are therefore the shear, S, which is accessed via the Maxwell and Reynolds stresses in the shearing 
wave, and the non-uniformity in the magnetic field, d z B, which is tapped through the radial electromotive force, discussed 
below. 



2.2 Dynamo loop and electromotive force 

The EMF is a fundamental quantity in dynamo theory (e.g. Mo ffatt 19781 . chap. 7). Whenever the velocity and magnetic 



fields are separated into mean and fluctuating parts, and the EMF is identified as the average of the cross product of the 
fluctuating velocity and magnetic fields, the curl of the EMF appears as a source term in the averaged induction equation. In 
our case, the EMF of the shearing wave, after azimuthal averaging, is 

= iRe(«* X 6). (12) 

If the nonlinear feedback of the shearing wave were considered, this EMF would generate a large-scale horizontal magnetic 
field B(z,t), which would evolve according to the azimuthally averaged induction equation 

d t B = B VC7 + V x£, (13) 
or, in components, 

d t B x = SB v -d z e v , (14) 

dtBy = d z £ x . (15) 

In the linear analysis that follows, we consider the waves to be of infinitesimal amplitude and the large-scale field B = B{z) e x 
to be constant in time. However, we are interested in identifying processes that could amplify this field or sustain it in the 
presence of non-zero resistivity. 

One possibility is that the radial EMF £ y directly amplifies the azimuthal field. This can be investigated by computing 
the correlation integral 

1=1 B x {-d z £ y )Az (16) 
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over the vertical extent of the system. Note that X is the rate of change of J \B 2 dz due to the radial EMF. If T > 0, the 
radial EMF acts to increase the energy of the large-scale azimuthal field, while the opposite is true if X < 0. Consistent with 
this, it can be seen from equation after an integration by parts that —I is also the rate at which the shearing wave derives 
energy from the large-scale field. We will see below that T is typically negative. 

A second possibility is that the azimuthal EMF £ x generates a radial field which, in combination with the shear, amplifies 
the azimuthal field. This can be checked using the correlation integral 

J= I B x (d z £ x )dz. (17) 



Assuming by convention that S > 0, we require J > for the second mechanism to work. It can be seen from equations (|14[) 
and (|15[) that, if J > 0, the B y generated by £ x is positively correlated with B x , so that when it is sheared it will amplify 
the existing B x . This method of obtaining a dynamo through a combination of the azimuthal EMF and the shear is related 
conceptually to the aQ. model of mean-field electrodynamics, but we do not think here of the EMF as arising from an a effect. 

In order to obtain a non-trivial EMF and non-zero values for the correlation integrals X and J , we must carry out the 
MRI calculation to a high degree of accuracy. Previous analytical treatments of the MRI have generally considered only a 
uniform magnetic field, or have relied on a local approximation or WKB method. These approaches are not accurate enough 
for our purposes, and result in a vanishing EMF because of the phase relationship between the velocity and magnetic field 
perturbations. Therefore we carry out a systematic asymptotic analysis of the MRI in the presence of a non-uniform magnetic 
field. 



2.3 Asymptotic analysis of localized solutions 

In the presence of a non-uniform magnetic field, the MRI will have preferred locations at which the growth rate of the instability, 
based on a local dispersion relation, is maximized. We therefore seek a solution that is localized in the neighbourhood of a 
favoured altitude z = z*. (More than one such position may exist, but we may consider the localized solutions independently.) 
To resolve the structure of the solution we introduce a stretched vertical coordinate £ such that 

z = z«+ eC, (18) 
where e 1 is a small dimensionless parameter that is used to organize an asymptotic expansion. We use a Taylor expansion 

°° n(«)™ 

B{z) = B, + BU( + \B'U 2 C + • • • = J] e n ^-±- (19) 
of the magnetic field, where the subscript * denotes evaluation at z — z*. The corresponding Alfven frequency is 

71 = 

The solution we are interested in has the form of a wavepacket, involving many wavelengths under a localizing envelope. 
It also grows exponentially in time to a first approximation, although this behaviour is modified on longer timescales. We 
introduce a slow time coordinate r such that 

t = e~\ (21) 

to describe the modulation of the exponential growth. 

The desired solution can be expressed as an asymptotic expansion in powers of e, having the form 

u x (z,t) = exp (yt + ie~ 3 k z z) [u x0 ((, r) + eu xl ((, r) + e 2 u x2 ((, r) H ], (22) 

and similarly for the other horizontal components u y , b x and b y , while 

u z (z,t) = eexp (7^ + ie~' i k z z) [u z0 ((,t) + ett z i(C, r) + e 2 u z2 ((,T) H ] , (23) 

and similarly for b z , and finally 

p(z,t) ^e 4 exp (jt + k z z) [ Po (<;,t) + e Pl ((,r) + e 2 P2 ((,T) + ■ ■ ■] (24) 

for the total pressure perturbation. This solution consists of a plane wave with wavevector (k x , —e~ 2 k x ST, e~ z k z ) and expo- 
nential growth rate 7, modulated by an envelope described by the functions u x o(C,,t), etc. The radial wavenumber is large 
when r = O(l) because the wave is then strongly sheared. We are interested in a regime in which the vertical wavenumber is 
larger still, here 0(e~ s ). Note that the various components of the solution are multiplied by different powers of e. The vertical 
components of the perturbations are smaller than the horizontal ones, while the vertical wavenumber is greater the radial one. 
This ordering allows the MRI to grow most rapidly, while maintaining consistency with the solenoidal nature of the velocity 
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and magnetic fields. As is usual in incompressible fluids, the pressure perturbation does whatever is required to maintain the 
solenoidal property of the velocity field. We assume formally that the horizontal perturbations are of order unity; since this 
is a linear theory, the overall scaling of the solution is arbitrary. 

We are interested in obtaining solutions in which the envelope is localized in the region where the coordinate £ is of order 
unity. This means that the physical scale of the localization in z is 0{e) while the wavelength in that direction is 0(e 3 ). The 
number of wavelengths contained within the envelope is then 0(e~ 2 ). Similarly, many e-foldings of the wave may occur within 
the timescale on which the envelope varies. 



2.4 Algebraic structure of the problem 



When the proposed asymptotic expansions are substituted into the linearized equations and the coefficients of each power of 
e are equated, we obtain a sequence of algebraic problems to be solved. At leading order, we obtain from equations © (|9} 



7 u *o = (20 — S)u v o + iw a *&xo, 
juyo = —2£lu x o + iu; a *bj,o, 

ju z0 = — ifc z po +i^a*bzO, 
-ybxo = Sbyo + iujatUxO, 

jbyO = lLJ^UyO, 

7& z0 = i^a*HzO, 

— ik x SrUyo + ik z u zQ = 0. 



(25) 
(26) 
(27) 
(28) 
(29) 
(30) 
(31) 



These equations are local and do not involve any derivatives with respect to £ because, as in the WKB theory, the vertical 
variation of the solution is determined at leading order by the rapid oscillation exp(ie _3 fc z z) under the modulatory envelope. 
The equations obtained at this and subsequent orders can be grouped systematically into a 'horizontal problem of order n\ 
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and a 'vertical problem of order n' 
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(33) 



where the forcing functions (A n ,...,G n ) depend on the solutions (u m ,b m ,p m ) at previous orders m < n, and on their 
derivatives with respect to £ and r. 

The first set of equations, n — 0, which are written out in full above, is homogeneous: Ao — Bo = Co = Do = Eq = 
Fo — Go = 0. In order to obtain a non-trivial solution, the determinant of the 4x4 matrix must vanish. This leads to the 
dispersion relation 

23( 7 , = 0, (34) 
where 

23( 7 , z) = ( 7 2 + ujI) 2 - 20S(7 2 + luI) + 4J7 V (35) 

is the usual dispersion relation for the MR0 in the limit that the vertical component of the wavevector dominates the other 
components. We assume that the system is hydrodynamically stable according to Rayleigh's criterion, i.e. Q(2£l — S) > 0. 
The MRI can then occur, with a positive real growth rate 7 > 0, if < to\ t < 2flS. Taking S > by convention, we will 
assume n > |S > (in a Keplerian disc, Q = (2/3)5). 

Subsequent sets of equations, n > 0, are inhomogeneous. Additional terms appear because of the Taylor expansion of 
B(z), the variation of the envelope of the wave with r and f, and other effects that are omitted at leading order. After some 
algebra, we find the following explicit expressions for the forcing functions: 



ik x p n —4 + ik x 



m— 1 



+ 



m— 1 



(m-1)! ' 



(36) 



1 Note, however, that in this dispersion relation w a is a function of z (see eq. 1201 1. 
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B 



E B (m) c .. 
— : by, n - m , (37) 
m! 



m—l 



ml ^ (m - 1)1 

m—l m — l 

Dn — Orby^n — 2 + lk x J Uy,n — mi (39) 



m — l 



-2 - d(p n -2 + ik x } — ■ b z , n -m, (40) 



m! 

m=l 



i^n = -cV&z,n-2 + ik x — ^ « Z ,n-m, (41) 



m! 



G„ = -ik x u X} „- 2 - d ( u z , n - 2 , (42) 

where it is understood that (it„, b n ,p n ) vanish for n < 0. 

Since the 4x4 matrix is singular when the dispersion relation is satisfied, the inhomogeneous horizontal problems can 
be solved only if 

2n-y 2 A n - 7 ( 7 2 + tol t )B n + 2Q,iiio^C n - ( 7 2 + loI, - 2£lS)\u^D n = 0. (43) 

We refer to this equation as the solvability condition of order n. It states that the forcing vector [A n B n C n D n ] T is 
orthogonal to the null eigenvector of the adjoint 4x4 matrix. 

Provided that this condition is met, the horizontal variables at 0(e n ) can be determined, up to the addition of a multiple 
of the null eigenvector of the 4x4 matrix. The solution of the horizontal problem of order n is 

_ 7-Bn + i^a*-Pn ~ ( 7 2 + uj,)Vn ( .., 

Uxn ~ 2th/ ' ( ' 

U yn = V n , (45) 

2ft 7 2 



, _ iuo^B n + 2QyC n + {2QS - ujl,)D n + i^ a *(2»S - 7 2 - ujj t )V n 

O xn — 7T7^ 2 ' l^O) 



Oyn — , 

7 

where the radial velocity V„(£,r), say, remains undetermined at this order. 
The solution of the vertical problem then follows as 

+ ik x SrV„ 
, k z F n + u!a. t G„ + iLO^kxSrVn 



-yk z 



(47) 

(48) 
(49) 



-ik z ~/E„ + ui & *k z F n + (7 2 + u>i»)G n + ( 7 2 + ul*)ik x STV n 
Pn = • (50) 

The determinant of the 3x3 matrix is 7 fc 2 and does not vanish. 

2.5 Development of the solution 

The first solvability condition to arise is that of order 1, 

4k x u; a 4-y 2 + lu'L - nS)B' t (V = 0. (51) 
which is equivalent to 

[a z O( 7 ,z»)]CVo = 0. (52) 

Physically, this relation states that the solution must be localized near an extremum of the local MRI growth rate and thereby 
allows us to determine z*. Since the solutions localized near the global maximum of the growth rate will dominate at later 
times, we will consider only them in the following. To find the possible locations, we note that the maximum possible MRI 
growth rate occurs for an optimal value of the Alfven frequency, given by cj 2 „ = S(Q ~ jS). Therefore, when this optimal 
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Figure 1. Example of localization for a sinusoidal magnetic profile. B(z) is plotted with a solid line and the local growth rate 7 
with a dashed line. The example here assumes S = 1, Q = 2/3 and k x = w. (i) (left): Weak magnetic field case. The preferred 
locations are the maxima of B 2 : z* = 0,0.5, 1. (ii) (centre): Strong magnetic field case. The preferred locations are the maxima of 7: 
z* = 0.13; 0.37; 0.63; 0.87. (iii) (right): Very strong magnetic field case. The preferred location is the minimum of B 2 : z t = 0.5. (iv) (not 
shown): Extremely strong magnetic field case: 7 = everywhere; no MRI appears in the system. 



value is found somewhere in the profile, the solution will be localized around this point. Otherwise, it will be localized near 
an extremum of B. According to these remarks, we find four non-trivial possibilities for the localization of the solution: 

(i) < maxuj < S(fl — jS), in which case the preferred location is at the maximum of B 2 {B'„ — solution in eq. I51|l : 

(ii) lu 2 = S(fl — jS) at some point, in which case that is the preferred location (j 2 = QS — lo 2 , solution in eg. 15 1 P : 

(iii) S(Q — jS) < mintj 2 < 2QS, in which case the preferred location is at the minimum of B 2 (B' t = solution in eg. 151}) : 

(iv) minima > 2Q.S, in which case the MRI does not occur. 

Examples of the various possible localizations for a sinusoidal profile are given in Fig. [1] 

The solvability condition of order n + 2 (for n > 0) involves the derivatives with respect to r of the horizontal variables 
of order n. It has the form 

AdrV n + (Br 2 + CC, 2 )V n = H n , (53) 
where 

A = d 7 V(-y,z*) 

= ^tf+Lul+n^Q-S)], (54) 
£ = S(7 2 +^,)V, (55) 

K, z 

c = Idlv^^z,) 

= 2k 2 x [(7 2 + 3ua» - Q,S)B* + (7 2 + w a 2 , - QS)B,B"] (56) 

are real constants, and H„ depends on the solutions V m at previous orders m < n. For a growing MRI solution, A and B are 
positive. 

We start with the case n = 0. Since Hq — 0, equation (I53f) is homogeneous and simply says that the local growth rate of 
the wave follows the local dispersion relation. Relative to the overall growth factor exp(7t), the envelope (described at leading 
order by the function Vb) decays slowly. The Bt 2 term reflects the effects on the local dispersion relation of a growing radial 
wavenumber, which reduces the local growth rate; the CC, 2 term corresponds to the vertical variation of the local growth rate 
around the extremum. The general solution is 

Vb(C,T) = W (C)£(C,T), (57) 
where 

*(f,r)=e*, (58) 

and W / o(C) is a function that depends on the initial conditions. Provided that C > 0, the wave becomes increasing vertically 
localized as time progresses. The condition C > means that the local MRI growth rate is maximized at z — z*, so that the 
solution becomes increasingly peaked at that location. 

Although the r 3 term in the exponential moderates the rapid exponential growth exp(e _2 7r) of the wave, it cannot be 
concluded from this analysis that the wave would ultimately decay. The asymptotic analysis does not apply in the limit of 
large r because it employs an ordering scheme in which the radial wavenumber is small compared to the vertical one. However, 
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it is true that, when the radial wavenumber has grown sufficiently to become comparable to or greater than the vertical one, 
the MRI will be weakened and ultimately suppressed. 

The asymptotic analysis also does not apply in the limit r — * 0, where the wave changes from leading to trailing and the 
condition \k y \ \k x \ does not hold. Therefore our solution has a typical 'intermediate asymptotic' character. Since we are 
interested in the total effect of one shearing wave during its whole lifetime, we focus only on trailing waves, as the maximum 
wave amplitude (and therefore maximum effect) is reached during this stage. Unless the wave is initialized when it is already 
strongly trailing, the function Wb(C) (and similar functions at higher orders) cannot be determined within this analysis, 
because the initial conditions cannot be applied within the scope of the solution. It may reasonably be assumed, however, that 
VKo(C) is a constant, because any vertical structure in the initial conditions is likely to be overwhelmed by the increasingly 
narrow Gaussian that develops according to the above solution. We therefore consider trailing waves and take Wo to be a 
constant, which may further be assumed to be real without loss of generality. When solving equation f|53p for n > 0, we do 
not add an arbitrary multiple of the complementary function E, since such terms depend on the initial conditions and can be 
absorbed into the definition of Wo- 

We have written a M ATHEMATICA script to solve the above systems of equations. It begins by expressing all the pertur- 
bation quantities (up to some specified truncation order) in terms of the variables V n and the forcing functions (A n , • ■ • , G n ) 
according to equations (144[) - (|50|) . It then evaluates the forcing functions (|36[) (|42l) . proceeding order by order so that eventu- 
ally all the perturbation quantities are expressed just in terms of the variables V n and their derivatives. Next, the solvability 
conditions (|53[) are determined by evaluating the functions H n . These equations are solved in turn (noting that E provides 
an integrating factor) to determine the explicit forms of V n (C, r) for n > 1. 



2.6 EMF and correlation integrals 

Our main interest is in calculating the horizontal EMF and its consequences for the large-scale magnetic field through the 
correlation integrals (|16p and (I17|) . Using the Taylor expansion (|19|l of B(z), we find 

J = <? exp(2 7 f)(2o + eTi + e 2 J 2 + ■••), (59) 
J = e 2 exp(2 7 t)(Jo + eji + e 2 J 2 + •••), (60) 
where 

f°° " D('"+ 1 )A"I 

ln= V~ ^£y,n- m d{, (61) 

/ L — ' m! 



' m=0 



, , g(m+l)/m 

Jn=~ > — r^—Sx.n-mdC, (62) 

/ ^ — ' m! 

and the horizontal EMF at order n (omitting a prefactor of eexp(27t)) is given by 

(£, T~) — ^Re I ^ ^ Uy m bz ,n — m ^zm^y,n — m J , (63) 

\m=0 / 

£yn(C,T) = |Re 2J "'m^,"-™ - y-*xrnbz,n~ m I . (64) 
\m=0 / 

The integration by parts used here, and the extension of the integrals to all values of £, are justified by the localized nature 
of the solution. 

An important property is that £ x o = £ y o = 0, which follows immediately from the phase relationships present in the 
leading-order solution. The Alfvenic couplings in the exponentially growing MRI shearing wave impose a phase shift of 7r/2 
between uo and bo- Therefore To = Jo = 0. 

In order to obtain a non-zero EMF, this phase relationship must be broken. The full expressions for H n and V n are very 
complicated, especially for larger n, and we will not write them here. The main interest is in their imaginary parts, which are 
partly responsible for introducing the required phase shift. We find 

Im(-ffi) =0 Im(Vi) = 0. (65) 

However, 

Im(# 2 ) = -y-W T 3 {E => Im(U 2 ) = -^-W r A C,E. (66) 
The expressions for the EMF at first and second orders are 



Localized MRI and dynamo 9 



£ x i = 0, (67) 
£ y , = -^0^\W o \ 2 r 2 E 2 , (68) 

^ 2 = ?£ §r |w,iVci?2 ' (69) 



(70) 



£ V 2 = ^fff| (T 2 + ^L)[2^. (7 2 + "a* - OS) - ( 7 2 + ^»)«S]|M/o| 2 r 2 C J B 2 + terms involving B'„ 

To proceed further we recall that there are two separate cases to consider, depending on whether B' t = or Bi 7^ 0. In 
the latter case, Bi 7^ 0, the dominant component of the EMF is £ y i and it has the correct symmetry to contribute to the 
correlation integral X\ . We then find 

/oo 
Bi£yi dC 
- 00 

klS^B? 2 2 (nA\^ 2 ( 2Br 3 \ 



2 7 fc 2 1 V2Cr7 " ^ 3.4 

The fact that Ii < means that the growing shearing wave acts to reduce the energy of the large-scale toroidal magnetic 
field. This effect can easily be understood as resulting from the vertical mixing of the non-uniform azimuthal field by the 
wave. It could be labelled a 'turbulent' resistivity. 

In the case B' t = 0, the leading contribution to T is 



Is = / B':C£ y2 d<Z 



= £^y+M^tf+>& - <*) - (V +«DmW fr (gi) 1/2 exp , (72) 

If the field is weak in the sense that tj^* is smaller than its optimal value S(Q — jS), then I3 < 0. In this case the wave still 
has a resistive effect, but a much weaker one. 

The leading contribution to the second correlation integral J is, in the general case, 

Js = - f (b':C£,2 + B',£ x3 \ dC, 

J —00 



(73) 

which simplifies to 

k 2 x S(B*B': + B?) ._„ l2 /nA\ 1/2 ( 2Bt s \ 

We note that this is proportional to — (B 2 )". As discussed in Section [2.21 we can identify a positive dynamo feedback when 
Jz > 0. When the field is weak and the wave is localized near a maximum of B 2 , we indeed have J3 > 0. On the other hand, 
when the field is very strong so that the wave is localized near a (non-zero) minimum of B 2 , we have < 0. When the wave 
is localized instead at a value of B that gives the optimal MRI growth rate, the outcome depends on the second derivative 
of B 2 at that point. A transition occurs between the case of a weak field, in which > 0, and the case of a strong field, in 
which Jz < 0. 

In Appendix [4] we give the details of some specific examples of localized shearing waves. 



3 COMPARISON WITH NUMERICAL SOLUTIONS OF SHEARING WAVES 
3.1 Comparison of high-fc z wavepackets 

We first compare the above linear results with a numerical calculation solving the linearized equations (f3)l- (|10[) . The numerical 
solution is computed assuming B = Bo cos(koz) with ko = 2tt/L z , k x = 2ir/L x , L x /L z = 2 and L z — 1, with a resolution 
of 512 Fourier modes in the vertical direction. During a numerical calculation, we evolve a single shearing wave, initialized 
at t = — 10 S^ 1 with a single Fourier mode k z = 125 ko. We then compute the electromotive force profiles at t = 50 S^ 1 and 
compare them with the predicted EMFs ( equations l67ffT0]i up to third ordei[f] for £ x and second order for £ y assuming e = 0.2 
and t — 2. To compare the results easily, and to remove the effect of the (arbitrary) amplitude of the initial perturbation, 



2 The third order for £ x is required to get the reversal described below. 
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Figure 2. Vertical profiles of the azimuthal and radial EMFs E x (left) and C y (right) for Bo = 0.1 (top) and Bo = 0.15 (bottom) at 
t = 50 S -1 with e = 0.2. As expected, the waves are localized near the maxima of the large-scale field at z = 0, 0.5 and f. 



the EMF profiles are normalized by the maximum of £ y . We plot in Fig. [5] examples of such comparisons in cases with weak 
fields (Bo < 0.2) and in Fig. [3] examples with strong fields (Bo > 0.2). 

These comparisons show that our analytical calculation is very close to the full linear computation. In the weak-field 
cases, we note that £ x and £ y have the same shape, as expected from (|69p and (|70p . This shape also naturally explains the 
resistive and dynamo effects quantified by integrals (|72[) and (|74[) . Moreover, the wavepackets are localized near the maxima 
of B at z = 0, 0.5 and 1, as expected in weak-field cases. In the strong-field cases, the wavepackets are localized near the 
maximum of 7 and, as expected, the resistive effect also becomes much stronger than the dynamo effect (\£ y \ ^> \£x\). We find 
a very good agreement between the two calculations, including a small asymmetry in £ x profile, resulting in the sign change 
of (for example, the double p eak at z ~ 0.35 h a s a la rger negative peak than a positive one). We note that the same kind 
of reversal has been observed by Lesur fc Qgilviei ( 20081 ) in the strong-field case, and was shown to be related to a possible 
dynamo cycle in accretion discs. 



3.2 Comparison of arbitrary waves 

To make a more general comparison, we c ompa re our analytical theory with the linear evolution of a random set of shearing 
waves, as computed by Lesur fc Qgilviei ( 20081 ) . We therefore use the same setup as in the previous case, except that we 
excite randomly all the vertical modes available. Very small viscosity and resistivity (y = r\ = 2 X 10~ 6 L^S) are added to 
damp short wavelengths. We then compute the EMF profiles at t — 30 S 1-1 for a set of 400 different initial conditions and 
compare the mean profiles with the analytical prediction for e = 0.2. Note that this choice of e is arbitrary because no single 
vertical wavenumber dominates in the numerical calculation. In the weak-field cases (Fig. [4|, the profiles are comparable 
to our asymptotic analysis, except near the peaks of £ x where some short-wavelength noise is observed. Nevertheless, these 
results show that the behaviour of the EMF in the weak-field cases is well described by our analysis, even in a more general 
situation where no assumption is made about the initial condition s , and small-fc z modes are allowed. Therefore, the positive 
dynamo effect and the resistive effect described by Lesur fc Ogilvie ( 20081 ) correspond closely to the integrals ([72} and (|74[) . In 
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Figure 3. Vertical profiles of the azimuthal and radial EMFs £ x (left) and £ y (right) for Bo = 0.3 (top) and Bo = 0.4 (bottom) at 
t = 50 S~ 1 with e = 0.2. The waves appear to be localized near the maxima of the growth rate, as in our analytical theory. 



the strong- field cases (Fig. [SJ, the comparisons are less satisfactory. We first note that the localization of the peaks does not 
correspond to the predictions of our analysis but occurs in regions with weaker field. This effect may occur because the local 
growth rate has a rather broad maximum with respect to z and the waves are not strongly localized, especially for smaller 
kg. We also find that the asymmetry in the double peak of £ x observed previously is enhanced in this case, showing that the 
reversal of the dynamo effect might be related to the correlations of waves with small k z . 



4 CONCLUSIONS 

In this paper, we have investigated the magnetorotational instability in the presence of an azimuthal magnetic field with a non- 
trivial vertical structure. We have shown first that, in the limit of perturbations of small vertical wavelength, the instability 
takes the form of wavepackets that are vertically localized near the local maxima of the MRI growth rate according to the 
standard dispersion relation. We have determined analytically the spatiotemporal evolution of these wavepackets. Computing 
the electromotive force that arises from the correlation of the velocity and magnetic field perturbations, we have investigated 
what would be the feedback of the MRI wavepackets on the large-scale field. We have found that the dominant process is 
a dissipative effect on the mean azimuthal field B x (z), which can be understood to a first approximation as a turbulent 
resistivity resulting from the mixing of the non-uniform azimuthal field. We have demonstrated that the MRI shearing waves 
can also have a feedback on the mean radial field B y (z), leading to the possibility of dynamo action when combined with 
the Keplerian shear. We have compared our analytical results with linear numerical calculations, showing that our analysis 
captures most of the physics involve d in these MRI shearin g waves. Finally, we have found that these results are consistent 
with the closure model presented by Lesur fc Qgilvie ( 20081 ) for nonlinear dynamo cycles in accretion discs. 



We emphasize that the findings presented here are only linear results that have many limitations. Hence, when some EMFs 
are referred to as providing a 'turbulent resistivity', one should understand that some energy is transferred from the large- 
scale field to the perturbation, but the total energy is conserved. When the system becomes nonlinear, these perturbations 
couple with themselves and the energy will eventually cascade towards small scales to be dissipated by molecular processes. 
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Figure 4. Vertical profiles of the azimuthal and radial EMFs E x (left) and C y (right) for Bo = 0.1 (top) and Bo = 0.15 (bottom) at 
t = 30 S 1- 1 with e = 0.2. The numerical profile is an average made from 400 shearing waves initialized at t = — 10 S~ 1 with random 
initial conditions. 



Therefore, the mechanism by which the energy is actually dissipated is intrinsically nonlinear, and is not covered by our work. 
Moreover, the amplitude of the quasilinear effects cannot be directly computed from our analysis, first because this amplitude 
will depend on the initial excitation of the waves, which is related to the strength of the underlying turbulent motions and is 
therefore unknown, and secondly because turbulent motions will modify the linear shearing waves as they grow and eventually 
destroy them. Since this destruction process will depend on the intensity of the underlying turbulence, the saturation level of 
the shearing waves is unpredictable, and may even be time-dependant. Therefore, our analysis can help us to understand the 
main effects due to the presence of the MRI in a turbulent flow, but it cannot predict the saturation values, which must be 
computed from nonlinear simulations. As an example, we cannot directly conclude from (|74[) that the dynamo effect should 
become stronger as we increase k x , because in a turbulent flow the nonlinear coupling will also become more efficient as we 
increase k. 

Despite these limitations, one can still deduce several important features from this linear analysis, the most important 
one being probably the dynamo feedback carried by £ x - This feedback allows us to close the loop given by the azimuthally 
averaged equations (|f4[) (|15[) . following a kind of af2 dynamo concept, although we do not consider the azimuthal EMF to 
arise from an a effect. Note however that since this feedback creates a large-scale radial field B y (z), this large-scale field should 
also be included in the linear analysis. We have chosen to neglect B y compared to B x , assuming that the dynamo feedback 
is smaller than the effect of the shear. Th is assumption is compatible with MRI turbulence in zero-net-flux simulations in 
which y/{Bi) > y/{Bl) jStone et al.lll99rj ). This feedback is therefore a natural explanation to the persistence of a magnetic 
field in zero-net-flux MRI simulations, but does not constitute on its own a self-sustaining mechanism, since one still needs 
an initial perturbation to excite non-axisymmetric waves. One may also note that our results are found in the ideal MHD 
limit, and therefore do not provide any information about the dependence of the feedback on the magnetic Prandtl number. 
To investigate this effect, we have tried introducing a small amount of resistivity and viscosity (< 10 -3 ) in our numerical 
calculations. It appears that the main results are not drastically modified by these effects, and no strong magnetic Prandtl 
number dependence is detectable from our linear analysis. This result may imply that the dynamo process itself does not 
depend strongly on the magnetic Prandtl number, although the way turbulence is maintained in the system does. This 
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conclusion is also compatible with simulations with a non-zero mean vertical flux IjLesur fc Longarettil 120071 ). in which no 
dynamo process is required to sustain the magnetic field but a strong magnetic Prandtl number dependence is still found. 

It has already be pointed out by iBrandenburg et al.l (119950 that MRI turbulence may act as an accretion disc dynamo. 
Using shearing box simulations with specific vertical boundary conditions, they have shown that large -scale azimuthal mag netic 
fields may be produced by MRI turbulence over long time scales. This idea has also been explored bv lRincon et al.l ([2007) with 
a stationary nonlinear so lution involving the M RI in a rotating plane Couette flow. More recently, we have found a dynamo 
cycle in MRI simulation i Lesur fc Qgilviel 20081 ) which motivated this linear analysis. Naturally, it would be tempting to say 
that all these numerical results are, at least partially, explained by our results. However, to confirm this conjecture, one has 
to reproduce all these numerical results, and check how the EMF profiles are correlated to the large-scale fields. If our result 
proves to be more general, it would be the first step toward a general closure model for global disc simulations, involving a 
self-consistent mechanism to generate magnetic fields. 
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APPENDIX A: ANALYTICAL EXAMPLES OF LOCALIZED WAVES 

The analytical analysis presented in this paper allows for a general vertical variation B(z) of the azimuthal field. Here we 
consider some simple examples. 

Let Q — (2/3)5', B(z) = Bo cos(2nz/L z ) and kx = ir/L z . This choice is motivated by our previous work (jLesur fc Qgilvie 
2008) and corresponds to a Keplerian shearing box in which the azimuthal field has a sinusoidal vertical dependence with 



wavelength L z , while the shearing wave of interest has an azimuthal wavelength of 2L Z . We then choose units such that 
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S = L z = 1. By also setting k z = 2-k, and recalling that the vertical wavenumber is e 3 fc z , we define the meaning of e to be 
the cube root of the ratio of the vertical wavelength of the shearing wave to that of the large-scale field. 

Al A weak- field example 

If B < \/5/(2\/37r) » 0.2055, then the wave is localized about z = z, — and other equivalent points at which ui 2 is 
maximized. This corresponds to case (i) in Fig. [T] 

As an example, let Bo = l/(2-\/37r). Then the basic growth rate is 7 = VE/6. We find 

/ r 3 tt 2 C 2 t\ 

E =^[-WE-^r)> (A1) 

= ^r\\W \ 2 E 2 , (A2) 
50 

|W ° |2 ■ : t „, x , )( — ,A1) 



£ y2 = ^^t\\W ?E 2 , (A3) 



J 3 = (V5r) exp ( ) . (A5) 



4 5 3 / 4 ^2lr7 V 36^5/ 

\W Q \ 2 , « . / r 3 

■ — ^=(Vor)exp -= 

4 5 3 / 4 V2^7 V 36^5 



A2 A strong-field example 

If B > %/5/(2 % /3tt), then the wave is preferentially localized about 
2tt V 12tt 2 B, 



2 = Z * = — arCCOS I m T 2D2 I ( A6 ) 



, 

2 



and other equivalent points at which tv 2 — 5/12, with a basic growth rate of 7 = 1/2. This corresponds to case (ii) in Fig. \T\ 
We find 



E = exp 



r>Ta (l^Bg -5)C a T 



48 16 



(A7) 



£ v i = -^(Wflg - 5) 1/2 t 2 |Wo| 2 B 2 , (A8) 
4V3 

£ x3 = — ^-=(127r 2 Bg - 5) 1/2 T-4r 2 + 12tt 2 (83 - 847r 2 B 2 )< 2 + 5tt 2 (12tt 2 B% - 5)CV 
6144^3 L 

+3tt 4 (425 - 10807r 2 Bg + 1447t 4 B 4 K 4 t] r 2 | W ! 2 B 2 , (A10) 
T x = --I=(12. 2 B 2 - 5) 1 / 2 r 3 / 2 |W | 2 exp (-^j , (All) 

Js = -^L(12^ 2 B ( 2 - 5)- 1/2 (5 - 6^ 2 B 2 )r 1/2 | W \ 2 exp ( ~ ) . (A12) 



.VST u ' v u/ " ' V - 1 

We see that J3 reverses sign and becomes negative for Bo > y/E/ (v6~7r) ~ 0.2906. 

These expressions give the contributions to 1 and J from a single localization point. The contributions from all equivalent 
localization points are of identical form. Although formally J3 diverges as Bo — + v / 5/(2v / 3tt), the description of the localization 
breaks down in this limit and the solution cannot be taken literally. 
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